Recall that the simple linear regression model is:
\[
y_i = \beta_0 + \beta_1 x_i + \epsilon_i \quad \text{for} ~ i=1,\ldots,n,
\]
where the error terms \(\epsilon_1,\ldots,\epsilon_n\) are mean zero, independent, and identically distributed random variables.
Let \(\hat{\beta}_1\) denote the parameter estimate of \(\beta_1\). The bias \(\hat{\beta}_1\) in estimating \(\beta_1\) is defined as the expected (i.e., mean) value of \(\hat{\beta}_1 - \beta_1\). A \(100(1-\alpha_1)\)% confidence interval for the slope \(\beta_1\) (assuming independent, normally distributed errors) is \(\hat{\beta}_1 \pm t_{\alpha_1/2,\text{df}} ~ \hat{\sigma}_{\hat{\beta}_1}\), where \(\text{df}\) is the residual (a.k.a., error) degrees of freedom, \(\hat{\sigma}_{\hat{\beta}_1}\) is the standard error of the slope estimate, and \(t_{\alpha_1/2,df}\) is the \(\alpha_1/2\) upper quantile of the \(t\)-distribution with \(\text{df}\) degrees of freedom. The coverage of a confidence interval estimator is the long-run relative frequency that the procedure generates intervals that contain the true value.
The regression coefficients \(\beta_0\) and \(\beta_1\) can be estimated using the method of least-squares. In R, if x and y are numeric vectors of the same length, the least-squares parameter estimates, standard errors, degrees of freedom, etc. are obtained using:
fm <- lm(y ~ x)
fm
summary(fm)
While theory provides results on the bias and coverage in simple linear regression models, the goal is to perform a simulation study to investigate:
- The bias in estimating the slope, and
- The coverage of the \(100(1-\alpha_1)\)% confidence interval estimator for the slope based on the usual normality assumption.
The following function conducts such a simulation study. The arguments to the function are:
regressor: Independent variable passed as a numeric vector of arbitrary length.
intercept: True intercept passed as a numeric vector of length one.
slope: True slope passed as a numeric vector of length one.
error_distribution: Distribution of the error term passed as a character vector of length one, equaling one of:
- “normal”: Error terms are independent and standard normally distributed.
- “chisq”: Error term are independently chi-squared distributed with 0.5 degrees of freedom and shifted to the left by 0.5 (so that the mean of the error term is zero). This has the effect of violating the normality assumption of the error terms \(\epsilon\)’s.
- “correlated”: Error term for observation i is normally distributed with mean \(0.2 z_i\) and variance 1, where \(z_i = (x_i-\bar{x})/s\) is the standardized value of regressor for observation i. This has the effect of making the error terms \(\epsilon\)’s correlated, thereby violating the independence assumption.
n_reps: Number of replications.
alpha1: Equals 1.0 minus the confidence coefficient for the confidence interval estimator of the slope. For example, for 95% confidence intervals on the slope, alpha1 = 0.05.
alpha2: Equals 1.0 minus the confidence coefficient for the confidence intervals of the bias and coverage probabilities. For example, for a 90% confidence interval on the bias, alpha2 = 0.10.
(x <- seq(-2, 2, length = 50))
[1] -2.00000000 -1.91836735 -1.83673469 -1.75510204 -1.67346939 -1.59183673 -1.51020408 -1.42857143
[9] -1.34693878 -1.26530612 -1.18367347 -1.10204082 -1.02040816 -0.93877551 -0.85714286 -0.77551020
[17] -0.69387755 -0.61224490 -0.53061224 -0.44897959 -0.36734694 -0.28571429 -0.20408163 -0.12244898
[25] -0.04081633 0.04081633 0.12244898 0.20408163 0.28571429 0.36734694 0.44897959 0.53061224
[33] 0.61224490 0.69387755 0.77551020 0.85714286 0.93877551 1.02040816 1.10204082 1.18367347
[41] 1.26530612 1.34693878 1.42857143 1.51020408 1.59183673 1.67346939 1.75510204 1.83673469
[49] 1.91836735 2.00000000
intercept <- 57
slope <- 112
true_y <- intercept + slope * x
plot(x, true_y, 'l')

error <- rnorm(length(x), sd = 10)
simulated_y <- true_y + error
plot(x, true_y, type = 'l')
points(x, simulated_y)

(lm_out <- lm(simulated_y ~ x))
Call:
lm(formula = simulated_y ~ x)
Coefficients:
(Intercept) x
59.32 111.89
(bias <- lm_out$coefficients[2] - slope)
x
-0.1086681
alpha <- .1
ci <- lm_out$coefficients[2] + c(-1, 1) * qt(1 - alpha, df = lm_out$df.residual) * summary(lm_out)$coefficients[2, "Std. Error"]
c(
lower = ci[1],
estimate = lm_out$coefficients[2],
upper = ci[2]
)
lower estimate.x upper
110.1580 111.8913 113.6246
library(tidyverse)
library(gganimate)
animation_data <- tibble(
x = rep(x, 10),
true_y = intercept + slope * x,
y = intercept + slope * x + rnorm(length(x), sd = 50),
simulation = rep(1:10, each = length(x) / 10)
)
animation_data %>%
ggplot(aes(x = x, y = y)) +
geom_point() +
geom_abline(slope = slope, intercept = intercept, col = "blue") +
geom_smooth(method = "lm", se = FALSE, col = "red") +
transition_states(states = simulation) +
theme_bw()

slope_simulation <- function(
regressor,
intercept,
slope,
error_distribution = "normal",
n_reps,
alpha1,
alpha2
) {
n <- length(regressor) # Number of observations
z <- (regressor - mean(regressor)) / sd(regressor) # Used for correlated errors
count_in <- 0
difference <- numeric(n_reps)
for (i in 1:n_reps) {
if (error_distribution == "normal") {
error <- rnorm(n)
} else if (error_distribution == "chisq") {
error <- rchisq(n, 0.5) - 0.5
} else if (error_distribution == "correlated") {
error <- rnorm(n, 0.2 * z, 1)
} else {
stop("Unknown error distribution")
}
response <- intercept + slope * regressor + error
fm <- lm(response ~ regressor)
# Extract slope estimate
fit_slope <- fm$coefficients[2]
# Calculate difference between slope estimate and slope
difference[i] <- fit_slope - slope
# Confidence interval on fit_slope
slope_conf <- fit_slope + c(-1, 1) * qt(1 - alpha1/2, fm$df.residual) * summary(fm)$coefficients[2, "Std. Error"]
count_in <- count_in + (slope_conf[1] < slope && slope < slope_conf[2])
}
# Calculate coverage
coverage <- count_in / n_reps
# Confidence interval on coverage
coverage_ci <- coverage + c(-1, 1) * qnorm(1 - alpha2/2) * sqrt((coverage * (1 - coverage))/n_reps)
# Bias estimate
bias <- mean(difference)
# Confidence interval on bias
bias_ci <- bias + c(-1, 1) * qnorm(1 - alpha2/2) * sd(difference) / sqrt(n_reps)
# Output
list(
coverage = c(
lower = coverage_ci[1],
estimate = coverage,
upper = coverage_ci[2]
),
bias = c(
lower = bias_ci[1],
estimate = bias,
upper = bias_ci[2]
)
)
}
For each iteration of the n_reps iterations of the simulation, randomly generate the response (i.e., dependent) vector using the supplied intercept, slope, regressor, and error term distribution. Fit the simple linear regression model and compute a \(100(1-\alpha_1)\)% confidence interval on the slope parameter. Record whether it contains the true slope. Also, record the difference between the slope estimate and its true value.
The proportion of times that the confidence interval contains the true slope is a point estimate of the its coverage. Theory says that the coverage should be \(1-\alpha_1\) when the usual regression assumptions are met. In addition to providing a point estimate of the coverage, you will provide a \(100(1-\alpha_2)\)% confidence interval on the coverage. (Use the normal approximation to the binomial, which is justified by the Central Limit Theorem since nreps is large.)
The average difference between the slope estimate and its true value is a point estimate of the bias. In addition to providing a point estimate of the bias, you will provide a \(100(1-\alpha_2)\)% confidence interval on the bias. (Again, the Central Limit Theorem is applicable.)
The function should return the following six elements:
- A point estimate of the bias in estimating the slope
- The lower bound of a \(100(1-\alpha_2)\)% confidence interval for the bias
- The upper bound of a \(100(1-\alpha_2)\)% confidence interval for the bias
- A point estimate of the coverage of the \(100(1-\alpha_1)\)% confidence interval estimator of the slope
- The lower bound of a \(100(1-\alpha_2)\)% confidence interval for the coverage
- The upper bound of a \(100(1-\alpha_2)\)% confidence interval for the coverage
Evaluate
Consider the following questions:
- Under what error terms is the estimator of the slope biased?
- Does the confidence interval estimator have the right coverage under the normally distributed error terms?
- What is the coverage under the chi-squared distributed error terms?
- What is the coverage under the correlated error terms?
- For the error terms in which the coverage is not correct, under what situation is it noticeable? When, if ever, does this coverage problem go away? What phenomenon makes it go away?
simulation <- function(
x,
intercept = 7,
slope = 0.5,
n_reps = 5000,
alpha1 = 0.10,
alpha2 = 0.05) {
for (error_dist in c("normal", "chisq", "correlated")) {
print(paste("Results for", error_dist, "errors:"))
results <- slope_simulation(x, intercept, slope, error_dist, n_reps, alpha1, alpha2)
print(paste0("Bias: ", results$bias[2], " (", results$bias[1], ", ", results$bias[3], ")"))
print(paste0("Coverage: ", results$coverage[2], " (", results$coverage[1], ", ", results$coverage[3], ")"))
}
}
simulation(x = seq(-2, 2, length = 50))
[1] "Results for normal errors:"
[1] "Bias: -0.000847634253262622 (-0.00418252079468957, 0.00248725228816433)"
[1] "Coverage: 0.9008 (0.89251422542325, 0.90908577457675)"
[1] "Results for chisq errors:"
[1] "Bias: 0.00181781368695764 (-0.00154552069960945, 0.00518114807352473)"
[1] "Coverage: 0.8946 (0.886088661926246, 0.903111338073754)"
[1] "Results for correlated errors:"
[1] "Bias: 0.167612312729427 (0.164251314366211, 0.170973311092642)"
[1] "Coverage: 0.5998 (0.586219840783989, 0.613380159216011)"
simulation(x = seq(-2, 2, length = 3))
[1] "Results for normal errors:"
[1] "Bias: 0.000966198337616723 (-0.00884107296014191, 0.0107734696353754)"
[1] "Coverage: 0.8968 (0.888367608976036, 0.905232391023964)"
[1] "Results for chisq errors:"
[1] "Bias: -0.00725001709345427 (-0.0167822938467396, 0.00228225965983101)"
[1] "Coverage: 0.9486 (0.942479509777187, 0.954720490222813)"
[1] "Results for correlated errors:"
[1] "Bias: 0.103372396296765 (0.0934498034634754, 0.113294989130054)"
[1] "Coverage: 0.8942 (0.885674433319529, 0.902725566680471)"
LS0tCnRpdGxlOiAiU2ltdWxhdGlvbiBTdHVkeSBmb3IgU2ltcGxlIExpbmVhciBSZWdyZXNzaW9uIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpgYGB7ciBzZXR1cCwgaW5jbHVkZT1GQUxTRX0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KGVjaG8gPSBUUlVFLCBjYWNoZSA9IFRSVUUsIG1lc3NhZ2UgPSBGQUxTRSkKYGBgCgpSZWNhbGwgdGhhdCB0aGUgc2ltcGxlIGxpbmVhciByZWdyZXNzaW9uIG1vZGVsIGlzOgoKJCQKeV9pID0gXGJldGFfMCArIFxiZXRhXzEgeF9pICsgXGVwc2lsb25faSBccXVhZCBcdGV4dHtmb3J9IH4gaT0xLFxsZG90cyxuLAokJAoKd2hlcmUgdGhlIGVycm9yIHRlcm1zICRcZXBzaWxvbl8xLFxsZG90cyxcZXBzaWxvbl9uJCBhcmUgbWVhbiB6ZXJvLCBpbmRlcGVuZGVudCwgCmFuZCBpZGVudGljYWxseSBkaXN0cmlidXRlZCByYW5kb20gdmFyaWFibGVzLgoKTGV0ICRcaGF0e1xiZXRhfV8xJCBkZW5vdGUgdGhlIHBhcmFtZXRlciBlc3RpbWF0ZSBvZiAkXGJldGFfMSQuIFRoZSBiaWFzCiRcaGF0e1xiZXRhfV8xJCBpbiBlc3RpbWF0aW5nICRcYmV0YV8xJCBpcyBkZWZpbmVkIGFzIHRoZSBleHBlY3RlZCAoaS5lLiwgbWVhbikKdmFsdWUgb2YgJFxoYXR7XGJldGF9XzEgLSBcYmV0YV8xJC4gQSAkMTAwKDEtXGFscGhhXzEpJCUgY29uZmlkZW5jZSBpbnRlcnZhbCBmb3IKdGhlIHNsb3BlICRcYmV0YV8xJCAoYXNzdW1pbmcgaW5kZXBlbmRlbnQsIG5vcm1hbGx5IGRpc3RyaWJ1dGVkIGVycm9ycykgaXMKJFxoYXR7XGJldGF9XzEgXHBtIHRfe1xhbHBoYV8xLzIsXHRleHR7ZGZ9fSB+IFxoYXR7XHNpZ21hfV97XGhhdHtcYmV0YX1fMX0kLAp3aGVyZSAkXHRleHR7ZGZ9JCBpcyB0aGUgcmVzaWR1YWwgKGEuay5hLiwgZXJyb3IpIGRlZ3JlZXMgb2YgZnJlZWRvbSwKJFxoYXR7XHNpZ21hfV97XGhhdHtcYmV0YX1fMX0kIGlzIHRoZSBzdGFuZGFyZCBlcnJvciBvZiB0aGUgc2xvcGUgZXN0aW1hdGUsIGFuZAokdF97XGFscGhhXzEvMixkZn0kIGlzIHRoZSAkXGFscGhhXzEvMiQgdXBwZXIgcXVhbnRpbGUgb2YgdGhlICR0JC1kaXN0cmlidXRpb24Kd2l0aCAkXHRleHR7ZGZ9JCBkZWdyZWVzIG9mIGZyZWVkb20uIFRoZSBjb3ZlcmFnZSBvZiBhIGNvbmZpZGVuY2UgaW50ZXJ2YWwKZXN0aW1hdG9yIGlzIHRoZSBsb25nLXJ1biByZWxhdGl2ZSBmcmVxdWVuY3kgdGhhdCB0aGUgcHJvY2VkdXJlIGdlbmVyYXRlcwppbnRlcnZhbHMgdGhhdCBjb250YWluIHRoZSB0cnVlIHZhbHVlLgoKVGhlIHJlZ3Jlc3Npb24gY29lZmZpY2llbnRzICRcYmV0YV8wJCBhbmQgJFxiZXRhXzEkIGNhbiBiZSBlc3RpbWF0ZWQgdXNpbmcgdGhlCm1ldGhvZCBvZiBsZWFzdC1zcXVhcmVzLiBJbiBSLCBpZiB4IGFuZCB5IGFyZSBudW1lcmljIHZlY3RvcnMgb2YgdGhlIHNhbWUKbGVuZ3RoLCB0aGUgbGVhc3Qtc3F1YXJlcyBwYXJhbWV0ZXIgZXN0aW1hdGVzLCBzdGFuZGFyZCBlcnJvcnMsIGRlZ3JlZXMgb2YKZnJlZWRvbSwgZXRjLiBhcmUgb2J0YWluZWQgdXNpbmc6CgpgYGB7ciwgZXZhbD1GQUxTRX0KZm0gPC0gbG0oeSB+IHgpCmZtCnN1bW1hcnkoZm0pCmBgYAoKV2hpbGUgdGhlb3J5IHByb3ZpZGVzIHJlc3VsdHMgb24gdGhlIGJpYXMgYW5kIGNvdmVyYWdlIGluIHNpbXBsZSBsaW5lYXIKcmVncmVzc2lvbiBtb2RlbHMsIHRoZSBnb2FsIGlzIHRvIHBlcmZvcm0gYSBzaW11bGF0aW9uIHN0dWR5IHRvIGludmVzdGlnYXRlOgoKKyBUaGUgYmlhcyBpbiBlc3RpbWF0aW5nIHRoZSBzbG9wZSwgYW5kCisgVGhlIGNvdmVyYWdlIG9mIHRoZSAkMTAwKDEtXGFscGhhXzEpJCUgY29uZmlkZW5jZSBpbnRlcnZhbCBlc3RpbWF0b3IgZm9yIHRoZQpzbG9wZSBiYXNlZCBvbiB0aGUgdXN1YWwgbm9ybWFsaXR5IGFzc3VtcHRpb24uCgpUaGUgZm9sbG93aW5nIGZ1bmN0aW9uIGNvbmR1Y3RzIHN1Y2ggYSBzaW11bGF0aW9uIHN0dWR5LiBUaGUgYXJndW1lbnRzIHRvIHRoZQpmdW5jdGlvbiBhcmU6CgorIGByZWdyZXNzb3JgOiBJbmRlcGVuZGVudCB2YXJpYWJsZSBwYXNzZWQgYXMgYSBudW1lcmljIHZlY3RvciBvZiBhcmJpdHJhcnkKbGVuZ3RoLgorIGBpbnRlcmNlcHRgOiBUcnVlIGludGVyY2VwdCBwYXNzZWQgYXMgYSBudW1lcmljIHZlY3RvciBvZiBsZW5ndGggb25lLgorIGBzbG9wZWA6IFRydWUgc2xvcGUgcGFzc2VkIGFzIGEgbnVtZXJpYyB2ZWN0b3Igb2YgbGVuZ3RoIG9uZS4KKyBgZXJyb3JfZGlzdHJpYnV0aW9uYDogRGlzdHJpYnV0aW9uIG9mIHRoZSBlcnJvciB0ZXJtIHBhc3NlZCBhcyBhIGNoYXJhY3Rlcgp2ZWN0b3Igb2YgbGVuZ3RoIG9uZSwgZXF1YWxpbmcgb25lIG9mOgogICAgKyAibm9ybWFsIjogRXJyb3IgdGVybXMgYXJlIGluZGVwZW5kZW50IGFuZCBzdGFuZGFyZCBub3JtYWxseSBkaXN0cmlidXRlZC4KICAgICsgImNoaXNxIjogRXJyb3IgdGVybSBhcmUgaW5kZXBlbmRlbnRseSBjaGktc3F1YXJlZCBkaXN0cmlidXRlZCB3aXRoIDAuNQogICAgZGVncmVlcyBvZiBmcmVlZG9tIGFuZCBzaGlmdGVkIHRvIHRoZSBsZWZ0IGJ5IDAuNSAoc28gdGhhdCB0aGUgbWVhbiBvZiB0aGUKICAgIGVycm9yIHRlcm0gaXMgemVybykuICBUaGlzIGhhcyB0aGUgZWZmZWN0IG9mIHZpb2xhdGluZyB0aGUgbm9ybWFsaXR5CiAgICBhc3N1bXB0aW9uIG9mIHRoZSBlcnJvciB0ZXJtcyAkXGVwc2lsb24kJ3MuCiAgICArICJjb3JyZWxhdGVkIjogRXJyb3IgdGVybSBmb3Igb2JzZXJ2YXRpb24gaSBpcyBub3JtYWxseSBkaXN0cmlidXRlZCB3aXRoCiAgICBtZWFuICQwLjIgel9pJCBhbmQgdmFyaWFuY2UgMSwgd2hlcmUgJHpfaSA9ICh4X2ktXGJhcnt4fSkvcyQgaXMgdGhlCiAgICBzdGFuZGFyZGl6ZWQgdmFsdWUgb2YgcmVncmVzc29yIGZvciBvYnNlcnZhdGlvbiBpLiAgVGhpcyBoYXMgdGhlIGVmZmVjdCBvZgogICAgbWFraW5nIHRoZSBlcnJvciB0ZXJtcyAkXGVwc2lsb24kJ3MgY29ycmVsYXRlZCwgdGhlcmVieSB2aW9sYXRpbmcgdGhlCiAgICBpbmRlcGVuZGVuY2UgYXNzdW1wdGlvbi4KKyBgbl9yZXBzYDogTnVtYmVyIG9mIHJlcGxpY2F0aW9ucy4KKyBgYWxwaGExYDogRXF1YWxzIDEuMCBtaW51cyB0aGUgY29uZmlkZW5jZSBjb2VmZmljaWVudCBmb3IgdGhlIGNvbmZpZGVuY2UKaW50ZXJ2YWwgZXN0aW1hdG9yIG9mIHRoZSBzbG9wZS4gRm9yIGV4YW1wbGUsIGZvciA5NSUgY29uZmlkZW5jZSBpbnRlcnZhbHMgb24KdGhlIHNsb3BlLCBhbHBoYTEgPSAwLjA1LgorIGBhbHBoYTJgOiBFcXVhbHMgMS4wIG1pbnVzIHRoZSBjb25maWRlbmNlIGNvZWZmaWNpZW50IGZvciB0aGUgY29uZmlkZW5jZQppbnRlcnZhbHMgb2YgdGhlIGJpYXMgYW5kIGNvdmVyYWdlIHByb2JhYmlsaXRpZXMuIEZvciBleGFtcGxlLCBmb3IgYSA5MCUKY29uZmlkZW5jZSBpbnRlcnZhbCBvbiB0aGUgYmlhcywgYWxwaGEyID0gMC4xMC4KCmBgYHtyfQooeCA8LSBzZXEoLTIsIDIsIGxlbmd0aCA9IDUwKSkKaW50ZXJjZXB0IDwtIDU3CnNsb3BlIDwtIDExMgp0cnVlX3kgPC0gaW50ZXJjZXB0ICsgc2xvcGUgKiB4CgpwbG90KHgsIHRydWVfeSwgJ2wnKQpgYGAKCgpgYGB7cn0KZXJyb3IgPC0gcm5vcm0obGVuZ3RoKHgpLCBzZCA9IDEwKQpzaW11bGF0ZWRfeSA8LSB0cnVlX3kgKyBlcnJvcgpgYGAKCgpgYGB7cn0KcGxvdCh4LCB0cnVlX3ksIHR5cGUgPSAnbCcpCnBvaW50cyh4LCBzaW11bGF0ZWRfeSkKYGBgCgpgYGB7cn0KKGxtX291dCA8LSBsbShzaW11bGF0ZWRfeSB+IHgpKQpgYGAKCmBgYHtyfQooYmlhcyA8LSBsbV9vdXQkY29lZmZpY2llbnRzWzJdIC0gc2xvcGUpCmBgYAoKYGBge3J9CmFscGhhIDwtIC4xCmNpIDwtIGxtX291dCRjb2VmZmljaWVudHNbMl0gKyBjKC0xLCAxKSAqIHF0KDEgLSBhbHBoYSwgZGYgPSBsbV9vdXQkZGYucmVzaWR1YWwpICogc3VtbWFyeShsbV9vdXQpJGNvZWZmaWNpZW50c1syLCAiU3RkLiBFcnJvciJdCmMoCiAgbG93ZXIgPSBjaVsxXSwKICBlc3RpbWF0ZSA9IGxtX291dCRjb2VmZmljaWVudHNbMl0sCiAgdXBwZXIgPSBjaVsyXQopCmBgYAoKYGBge3IsIGV2YWwgPSBGQUxTRX0KbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkoZ2dhbmltYXRlKQoKYW5pbWF0aW9uX2RhdGEgPC0gdGliYmxlKAogIHggPSByZXAoeCwgMTApLAogIHRydWVfeSA9IGludGVyY2VwdCArIHNsb3BlICogeCwKICB5ID0gaW50ZXJjZXB0ICsgc2xvcGUgKiB4ICsgcm5vcm0obGVuZ3RoKHgpLCBzZCA9IDUwKSwKICBzaW11bGF0aW9uID0gcmVwKDE6MTAsIGVhY2ggPSBsZW5ndGgoeCkgLyAxMCkKKSAKCmFuaW1hdGlvbl9kYXRhICU+JSAKICBnZ3Bsb3QoYWVzKHggPSB4LCB5ID0geSkpICsKICBnZW9tX3BvaW50KCkgKwogIGdlb21fYWJsaW5lKHNsb3BlID0gc2xvcGUsIGludGVyY2VwdCA9IGludGVyY2VwdCwgY29sID0gImJsdWUiKSArCiAgZ2VvbV9zbW9vdGgobWV0aG9kID0gImxtIiwgc2UgPSBGQUxTRSwgY29sID0gInJlZCIpICsKICB0cmFuc2l0aW9uX3N0YXRlcyhzdGF0ZXMgPSBzaW11bGF0aW9uKSArCiAgdGhlbWVfYncoKQpgYGAKCiFbXShyZWdyZXNzaW9uX2FuaW1hdGlvbi5naWYpCgpgYGB7cn0Kc2xvcGVfc2ltdWxhdGlvbiA8LSBmdW5jdGlvbigKICByZWdyZXNzb3IsCiAgaW50ZXJjZXB0LAogIHNsb3BlLAogIGVycm9yX2Rpc3RyaWJ1dGlvbiA9ICJub3JtYWwiLAogIG5fcmVwcywKICBhbHBoYTEsCiAgYWxwaGEyCikgewogIG4gPC0gbGVuZ3RoKHJlZ3Jlc3NvcikgIyBOdW1iZXIgb2Ygb2JzZXJ2YXRpb25zCiAgeiA8LSAocmVncmVzc29yIC0gbWVhbihyZWdyZXNzb3IpKSAvIHNkKHJlZ3Jlc3NvcikgIyBVc2VkIGZvciBjb3JyZWxhdGVkIGVycm9ycwogIGNvdW50X2luIDwtIDAKICBkaWZmZXJlbmNlIDwtIG51bWVyaWMobl9yZXBzKQogIGZvciAoaSBpbiAxOm5fcmVwcykgewogICAgaWYgKGVycm9yX2Rpc3RyaWJ1dGlvbiA9PSAibm9ybWFsIikgewogICAgICBlcnJvciA8LSBybm9ybShuKQogICAgfSBlbHNlIGlmIChlcnJvcl9kaXN0cmlidXRpb24gPT0gImNoaXNxIikgewogICAgICBlcnJvciA8LSByY2hpc3EobiwgMC41KSAtIDAuNQogICAgfSBlbHNlIGlmIChlcnJvcl9kaXN0cmlidXRpb24gPT0gImNvcnJlbGF0ZWQiKSB7CiAgICAgIGVycm9yIDwtIHJub3JtKG4sIDAuMiAqIHosIDEpCiAgICB9IGVsc2UgewogICAgICBzdG9wKCJVbmtub3duIGVycm9yIGRpc3RyaWJ1dGlvbiIpCiAgICB9CiAgICAKICAgIHJlc3BvbnNlIDwtIGludGVyY2VwdCArIHNsb3BlICogcmVncmVzc29yICsgZXJyb3IKICAgIGZtIDwtIGxtKHJlc3BvbnNlIH4gcmVncmVzc29yKQogICAgCiAgICAjIEV4dHJhY3Qgc2xvcGUgZXN0aW1hdGUKICAgIGZpdF9zbG9wZSA8LSBmbSRjb2VmZmljaWVudHNbMl0KICAgIAogICAgIyBDYWxjdWxhdGUgZGlmZmVyZW5jZSBiZXR3ZWVuIHNsb3BlIGVzdGltYXRlIGFuZCBzbG9wZQogICAgZGlmZmVyZW5jZVtpXSA8LSBmaXRfc2xvcGUgLSBzbG9wZQogICAgCiAgICAjIENvbmZpZGVuY2UgaW50ZXJ2YWwgb24gZml0X3Nsb3BlCiAgICBzbG9wZV9jb25mIDwtIGZpdF9zbG9wZSArIGMoLTEsIDEpICogcXQoMSAtIGFscGhhMS8yLCBmbSRkZi5yZXNpZHVhbCkgKiBzdW1tYXJ5KGZtKSRjb2VmZmljaWVudHNbMiwgIlN0ZC4gRXJyb3IiXQogICAgY291bnRfaW4gPC0gY291bnRfaW4gKyAoc2xvcGVfY29uZlsxXSA8IHNsb3BlICYmIHNsb3BlIDwgc2xvcGVfY29uZlsyXSkKICB9CiAgCiAgIyBDYWxjdWxhdGUgY292ZXJhZ2UKICBjb3ZlcmFnZSA8LSBjb3VudF9pbiAvIG5fcmVwcwogIAogICMgQ29uZmlkZW5jZSBpbnRlcnZhbCBvbiBjb3ZlcmFnZQogIGNvdmVyYWdlX2NpIDwtIGNvdmVyYWdlICsgYygtMSwgMSkgKiBxbm9ybSgxIC0gYWxwaGEyLzIpICogc3FydCgoY292ZXJhZ2UgICogKDEgLSBjb3ZlcmFnZSkpL25fcmVwcykKICAKICAjIEJpYXMgZXN0aW1hdGUKICBiaWFzIDwtIG1lYW4oZGlmZmVyZW5jZSkKICAKICAjIENvbmZpZGVuY2UgaW50ZXJ2YWwgb24gYmlhcwogIGJpYXNfY2kgPC0gYmlhcyArIGMoLTEsIDEpICogcW5vcm0oMSAtIGFscGhhMi8yKSAqIHNkKGRpZmZlcmVuY2UpIC8gc3FydChuX3JlcHMpCiAgCiAgIyBPdXRwdXQKICBsaXN0KAogICAgY292ZXJhZ2UgPSBjKAogICAgICBsb3dlciA9IGNvdmVyYWdlX2NpWzFdLAogICAgICBlc3RpbWF0ZSA9IGNvdmVyYWdlLAogICAgICB1cHBlciA9IGNvdmVyYWdlX2NpWzJdCiAgICApLAogICAgYmlhcyA9IGMoCiAgICAgIGxvd2VyID0gYmlhc19jaVsxXSwKICAgICAgZXN0aW1hdGUgPSBiaWFzLAogICAgICB1cHBlciA9IGJpYXNfY2lbMl0KICAgICkKICApCn0KYGBgCgpGb3IgZWFjaCBpdGVyYXRpb24gb2YgdGhlIGBuX3JlcHNgIGl0ZXJhdGlvbnMgb2YgdGhlIHNpbXVsYXRpb24sIHJhbmRvbWx5CmdlbmVyYXRlIHRoZSByZXNwb25zZSAoaS5lLiwgZGVwZW5kZW50KSB2ZWN0b3IgdXNpbmcgdGhlIHN1cHBsaWVkIGludGVyY2VwdCwKc2xvcGUsIHJlZ3Jlc3NvciwgYW5kIGVycm9yIHRlcm0gZGlzdHJpYnV0aW9uLiBGaXQgdGhlIHNpbXBsZSBsaW5lYXIgcmVncmVzc2lvbgptb2RlbCBhbmQgY29tcHV0ZSBhICQxMDAoMS1cYWxwaGFfMSkkJSBjb25maWRlbmNlIGludGVydmFsIG9uIHRoZSBzbG9wZQpwYXJhbWV0ZXIuIFJlY29yZCB3aGV0aGVyIGl0IGNvbnRhaW5zIHRoZSB0cnVlIHNsb3BlLiBBbHNvLCByZWNvcmQgdGhlCmRpZmZlcmVuY2UgYmV0d2VlbiB0aGUgc2xvcGUgZXN0aW1hdGUgYW5kIGl0cyB0cnVlIHZhbHVlLgoKVGhlIHByb3BvcnRpb24gb2YgdGltZXMgdGhhdCB0aGUgY29uZmlkZW5jZSBpbnRlcnZhbCBjb250YWlucyB0aGUgdHJ1ZSBzbG9wZSBpcwphIHBvaW50IGVzdGltYXRlIG9mIHRoZSBpdHMgY292ZXJhZ2UuIFRoZW9yeSBzYXlzIHRoYXQgdGhlIGNvdmVyYWdlIHNob3VsZCBiZQokMS1cYWxwaGFfMSQgd2hlbiB0aGUgdXN1YWwgcmVncmVzc2lvbiBhc3N1bXB0aW9ucyBhcmUgbWV0LiBJbiBhZGRpdGlvbiB0bwpwcm92aWRpbmcgYSBwb2ludCBlc3RpbWF0ZSBvZiB0aGUgY292ZXJhZ2UsIHlvdSB3aWxsIHByb3ZpZGUgYQokMTAwKDEtXGFscGhhXzIpJCUgY29uZmlkZW5jZSBpbnRlcnZhbCBvbiB0aGUgY292ZXJhZ2UuIChVc2UgdGhlIG5vcm1hbAphcHByb3hpbWF0aW9uIHRvIHRoZSBiaW5vbWlhbCwgd2hpY2ggaXMganVzdGlmaWVkIGJ5IHRoZSBDZW50cmFsIExpbWl0IFRoZW9yZW0Kc2luY2UgbnJlcHMgaXMgbGFyZ2UuKQoKVGhlIGF2ZXJhZ2UgZGlmZmVyZW5jZSBiZXR3ZWVuIHRoZSBzbG9wZSBlc3RpbWF0ZSBhbmQgaXRzIHRydWUgdmFsdWUgaXMgYSBwb2ludAplc3RpbWF0ZSBvZiB0aGUgYmlhcy4gSW4gYWRkaXRpb24gdG8gcHJvdmlkaW5nIGEgcG9pbnQgZXN0aW1hdGUgb2YgdGhlIGJpYXMsIHlvdQp3aWxsIHByb3ZpZGUgYSAkMTAwKDEtXGFscGhhXzIpJCUgY29uZmlkZW5jZSBpbnRlcnZhbCBvbiB0aGUgYmlhcy4gKEFnYWluLCB0aGUKQ2VudHJhbCBMaW1pdCBUaGVvcmVtIGlzIGFwcGxpY2FibGUuKQoKVGhlIGZ1bmN0aW9uIHNob3VsZCByZXR1cm4gdGhlIGZvbGxvd2luZyBzaXggZWxlbWVudHM6CgoxLiBBIHBvaW50IGVzdGltYXRlIG9mIHRoZSBiaWFzIGluIGVzdGltYXRpbmcgdGhlIHNsb3BlCjIuIFRoZSBsb3dlciBib3VuZCBvZiBhICQxMDAoMS1cYWxwaGFfMikkJSBjb25maWRlbmNlIGludGVydmFsIGZvciB0aGUgYmlhcwozLiBUaGUgdXBwZXIgYm91bmQgb2YgYSAkMTAwKDEtXGFscGhhXzIpJCUgY29uZmlkZW5jZSBpbnRlcnZhbCBmb3IgdGhlIGJpYXMKNC4gQSBwb2ludCBlc3RpbWF0ZSBvZiB0aGUgY292ZXJhZ2Ugb2YgdGhlICQxMDAoMS1cYWxwaGFfMSkkJSBjb25maWRlbmNlCmludGVydmFsIGVzdGltYXRvciBvZiB0aGUgc2xvcGUKNS4gVGhlIGxvd2VyIGJvdW5kIG9mIGEgJDEwMCgxLVxhbHBoYV8yKSQlIGNvbmZpZGVuY2UgaW50ZXJ2YWwgZm9yIHRoZSBjb3ZlcmFnZQo2LiBUaGUgdXBwZXIgYm91bmQgb2YgYSAkMTAwKDEtXGFscGhhXzIpJCUgY29uZmlkZW5jZSBpbnRlcnZhbCBmb3IgdGhlIGNvdmVyYWdlCgojIyBFdmFsdWF0ZQpDb25zaWRlciB0aGUgZm9sbG93aW5nIHF1ZXN0aW9uczoKCisgVW5kZXIgd2hhdCBlcnJvciB0ZXJtcyBpcyB0aGUgZXN0aW1hdG9yIG9mIHRoZSBzbG9wZSBiaWFzZWQ/CisgRG9lcyB0aGUgY29uZmlkZW5jZSBpbnRlcnZhbCBlc3RpbWF0b3IgaGF2ZSB0aGUgcmlnaHQgY292ZXJhZ2UgdW5kZXIgdGhlCm5vcm1hbGx5IGRpc3RyaWJ1dGVkIGVycm9yIHRlcm1zPworIFdoYXQgaXMgdGhlIGNvdmVyYWdlIHVuZGVyIHRoZSBjaGktc3F1YXJlZCBkaXN0cmlidXRlZCBlcnJvciB0ZXJtcz8KKyBXaGF0IGlzIHRoZSBjb3ZlcmFnZSB1bmRlciB0aGUgY29ycmVsYXRlZCBlcnJvciB0ZXJtcz8KKyBGb3IgdGhlIGVycm9yIHRlcm1zIGluIHdoaWNoIHRoZSBjb3ZlcmFnZSBpcyBub3QgY29ycmVjdCwgdW5kZXIgd2hhdCBzaXR1YXRpb24KaXMgaXQgbm90aWNlYWJsZT8gV2hlbiwgaWYgZXZlciwgZG9lcyB0aGlzIGNvdmVyYWdlIHByb2JsZW0gZ28gYXdheT8gV2hhdApwaGVub21lbm9uIG1ha2VzIGl0IGdvIGF3YXk/CgpgYGB7cn0Kc2ltdWxhdGlvbiA8LSBmdW5jdGlvbigKICB4LCAKICBpbnRlcmNlcHQgPSA3LCAKICBzbG9wZSA9IDAuNSwgCiAgbl9yZXBzID0gNTAwMCwgCiAgYWxwaGExID0gMC4xMCwgCiAgYWxwaGEyID0gMC4wNSkgewogIGZvciAoZXJyb3JfZGlzdCBpbiBjKCJub3JtYWwiLCAiY2hpc3EiLCAiY29ycmVsYXRlZCIpKSB7CiAgICBwcmludChwYXN0ZSgiUmVzdWx0cyBmb3IiLCBlcnJvcl9kaXN0LCAiZXJyb3JzOiIpKQogICAgcmVzdWx0cyA8LSBzbG9wZV9zaW11bGF0aW9uKHgsIGludGVyY2VwdCwgc2xvcGUsIGVycm9yX2Rpc3QsIG5fcmVwcywgYWxwaGExLCBhbHBoYTIpCiAgICBwcmludChwYXN0ZTAoIkJpYXM6ICAgICAiLCByZXN1bHRzJGJpYXNbMl0sICIgKCIsIHJlc3VsdHMkYmlhc1sxXSwgIiwgIiwgcmVzdWx0cyRiaWFzWzNdLCAiKSIpKQogICAgcHJpbnQocGFzdGUwKCJDb3ZlcmFnZTogIiwgcmVzdWx0cyRjb3ZlcmFnZVsyXSwgIiAoIiwgcmVzdWx0cyRjb3ZlcmFnZVsxXSwgIiwgIiwgcmVzdWx0cyRjb3ZlcmFnZVszXSwgIikiKSkKICB9Cn0KCnNpbXVsYXRpb24oeCA9IHNlcSgtMiwgMiwgbGVuZ3RoID0gNTApKQpgYGAKCmBgYHtyfQpzaW11bGF0aW9uKHggPSBzZXEoLTIsIDIsIGxlbmd0aCA9IDMpKQpgYGAKCg==